House Prices - Advanced Regression Techniques¶
In [1]:
# Basic Libraries to import
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import math
from scipy.stats import norm, skew
import os
from quickda.explore_data import *
import warnings
warnings.filterwarnings("ignore")
%matplotlib inline
In [2]:
pd.pandas.set_option('display.max_columns',None)
pd.pandas.set_option('display.max_rows',90)
pd.options.display.float_format = '{:.5f}'.format
# Setting seaborn style
sns.set_palette('Reds')
sns.set_style('darkgrid')
In [3]:
!pip -q install kaggle
from kaggle.api.kaggle_api_extended import KaggleApi
api = KaggleApi()
api.authenticate()
In [4]:
api.competitions_list(category='gettingStarted')
Out[4]:
In [5]:
api.competition_download_files('house-prices-advanced-regression-techniques')
In [6]:
from zipfile import ZipFile
zf = ZipFile('house-prices-advanced-regression-techniques.zip')
zf.extractall('data\\')
zf.close()
In [7]:
housePrices = pd.read_csv('./data/train.csv')
train_index = housePrices['Id']
target = housePrices['SalePrice']
target_col = 'SalePrice'
housePrices.head()
Out[7]:
In [8]:
housePrices.shape
Out[8]:
In [9]:
housePrices.info()
In [10]:
test_df = pd.read_csv('./data/test.csv')
test_index = test_df['Id']
test_df.head()
Out[10]:
In [11]:
housePrices.shape, test_df.shape
Out[11]:
In [12]:
from IPython.display import display
with pd.option_context('display.max_rows', 100, 'display.max_columns', 100):
display(explore(housePrices, method="summarize"))
Missing value analysis¶
In [13]:
df = pd.concat((housePrices, test_df))
indexs = df['Id']
df.drop(['SalePrice', 'Id'],axis = 1, inplace = True)
df.shape, housePrices.shape
Out[13]:
In [14]:
df.head()
Out[14]:
In [15]:
percent_missing = round(df.isnull().sum() * 100 / len(df),3)
missing_value_df = pd.DataFrame({'column_name': df.columns,
'percent_missing': percent_missing})
missing_value_df[missing_value_df['percent_missing']>0].sort_values(by='percent_missing', ascending= False )
Out[15]:
Let's analyse whether the feature has:¶
- Missing value by random
- Missing value by chance
From the metadata and data description, it looks like the houses without Pool, MiscFeatures, Alley access, Fence and Fireplace are are marked as missing. Hence replace those values with 'NotAvailable'¶
In [16]:
colToReplace = ['PoolQC', 'MiscFeature', 'Alley', 'Fence', 'FireplaceQu']
for i in colToReplace:
df[i] = df[i].fillna('NotAvailable')
In [17]:
df['PoolQC'].value_counts()
Out[17]:
In [18]:
round(df[colToReplace].isnull().sum() * 100 / len(df),3)
Out[18]:
Garage¶
Houses with no Garage are missing. Hence replacing missing with 'NoGarage'
In [19]:
garageCols = ['GarageType', 'GarageFinish', 'GarageQual', 'GarageCond','GarageYrBlt', 'GarageArea', 'GarageCars']
for i in garageCols[0:4]:
df[i] = df[i].fillna('NoGarage')
In [20]:
round(df[garageCols].isnull().sum() * 100 / len(df),3)
Out[20]:
Replacing 'GarageYrBlt', 'GarageArea', 'GarageCars' with 0, since there is no garage.¶
In [21]:
for i in garageCols[4:]:
df[i] = df[i].fillna(0)
Basement¶
Houses with no basement are missing. Hence replacing missing with 'NoBasement'
In [22]:
basementCols = ['BsmtQual', 'BsmtCond', 'BsmtExposure', 'BsmtFinType1','BsmtFinType2']
round(df[basementCols].isnull().sum() * 100 / len(df),3)
Out[22]:
In [23]:
for i in basementCols:
df[i] = df[i].fillna('NoBasement')
round(df[basementCols].isnull().sum() * 100 / len(df),3)
Out[23]:
In [24]:
for i in ['BsmtFinSF1','BsmtFinSF2', 'BsmtUnfSF', 'TotalBsmtSF', 'BsmtFullBath','BsmtHalfBath']:
df[i] = df[i].fillna(0)
Masonry veneer¶
In [25]:
df['MasVnrType'] = df['MasVnrType'].fillna('None')
df['MasVnrArea'] = df['MasVnrArea'].fillna(0)
Functional¶
In [26]:
df['Functional'] = df['Functional'].fillna('Typ')
Lot Frontage¶
In [27]:
df['Neighborhood'].value_counts()
Out[27]:
Usually houses in same neighborhood look similar and therefore imputing missing values by median of lot frontage grouped by neighborhood¶
In [28]:
df['LotFrontage'] = df.groupby('Neighborhood')['LotFrontage'].transform(lambda x: x.fillna(x.median()))
In [29]:
percent_missing = round(df.isnull().sum() * 100 / len(df),3)
missing_value_df = pd.DataFrame({'column_name': df.columns,
'percent_missing': percent_missing})
missing_value_df[missing_value_df['percent_missing']>0].sort_values(by='percent_missing', ascending= False )
Out[29]:
In [30]:
cols = ['MSZoning', 'Utilities', 'Exterior1st', 'Exterior2nd', 'KitchenQual', 'SaleType', 'Electrical']
In [31]:
for i in cols:
df[i] = df[i].fillna(df[i].mode()[0])
In [32]:
plt.figure(figsize=(7,5))
sns.heatmap(df.isnull(), cbar=True, cmap= 'YlGnBu')
Out[32]:
Feature Engineering¶
In [33]:
categorical_cols = df.dtypes[ df.dtypes == 'object' ].index.tolist()
numerical_cols = df.dtypes[ df.dtypes != 'object' ].index.tolist()
print("Total Number of Features:", df.shape[1])
print("Total Number of Categorical Features:", len(categorical_cols))
print("Total Number of Numerical Features:", len(numerical_cols))
In [34]:
# Custom Encoding
MSSubclass_map = {'20':1,'30':2,'40':3,'45':4,'50':5,'60':6,'70':7,'75':8,'80':9,
'85':10, '90':11,'120':12,'150':13,'160':14,'180':15,'190':16}
LotShape = {'Reg':3,'IR1':2,'IR2':1,'IR3':0}
utilities_map = {'AllPub':4, 'NoSewr':3, 'NoSeWa':2, 'ELO':1}
quality_map = {'Ex':5,'Gd':4,'TA':3,'Fa':2,'Po':1, 'NoBasement':0, 'NotAvailable':0, 'NoGarage':0}
basement_exp = {'NoBasement':0,'No':1,'Mn':2,'Av':3,'Gd':4}
bsmtfin = {'NoBasement':0,'Unf':1,'LwQ':2, 'Rec':3,'BLQ':4,'ALQ':5, 'GLQ':6}
centralAir = {'N':0, 'Y':1}
Functional_map = {'Sal':0,'Sev':1,'Maj2':2,'Maj1':3,'Mod':4,'Min2':5,'Min1':6,'Typ':7}
GarageFinish = {'NoGarage':0,'Unf':1,'RFn':2,'Fin':3}
fence_qc = {'GdPrv':4, 'MnPrv':3, 'GdWo':2, 'MnWw':1, 'NotAvailable':0}
In [35]:
q = ['ExterQual', 'ExterCond', 'BsmtQual', 'BsmtCond',
'KitchenQual', 'FireplaceQu', 'GarageQual', 'GarageCond', 'PoolQC', 'HeatingQC']
for i in q:
df[i] = df[i].replace(quality_map)
In [36]:
df['MSSubClass'] = df['MSSubClass'].apply(str).replace(MSSubclass_map)
df['LotShape'] = df['LotShape'].replace(LotShape)
df['Utilities'] = df['Utilities'].replace(utilities_map)
df['BsmtExposure'] = df['BsmtExposure'].replace(basement_exp)
df['BsmtFinType1'] = df['BsmtFinType1'].replace(bsmtfin)
df['BsmtFinType2'] = df['BsmtFinType2'].replace(bsmtfin)
df['CentralAir'] = df['CentralAir'].replace(centralAir)
df['Functional'] = df['Functional'].replace(Functional_map)
df['GarageFinish'] = df['GarageFinish'].replace(GarageFinish)
df['Fence'] = df['Fence'].replace(fence_qc)
In [37]:
#Creating a feature -> If house is remodelled then 1 else 0
df['Remodel'] = np.where(df['YearRemodAdd']==df['YearBuilt'], 1, 0)
df.drop('YearRemodAdd', axis=1, inplace=True)
In [38]:
#Transforming variable to age
df['HouseAge'] = 2020 - df['YearBuilt']
df.drop('YearBuilt', axis=1, inplace=True)
In [39]:
#Transforming variable to age
df['GarageYrBlt'] = 2020 - df['GarageYrBlt']
df.drop('GarageYrBlt', axis=1, inplace=True)
In [40]:
df['TotalSF']=df['TotalBsmtSF'] + df['1stFlrSF'] + df['2ndFlrSF']
df['Total_Bathrooms'] = (df['FullBath'] + (0.5 * df['HalfBath'])+df['BsmtFullBath']+(0.5 * df['BsmtHalfBath']))
df['Total_porch_sf'] = (df['OpenPorchSF'] + df['3SsnPorch'] +df['EnclosedPorch'] +
df['ScreenPorch'] + df['WoodDeckSF'])
In [41]:
# Import library for VIF
from statsmodels.stats.outliers_influence import variance_inflation_factor
def calc_vif(X):
# Calculating VIF
vif = pd.DataFrame()
vif["variables"] = X.columns
vif["VIF"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
return(vif)
def annote_graph(graph):
for p in graph.patches:
height = p.get_height()
graph.text(p.get_x()+p.get_width()/2., height + 0.1,height ,ha="center")
In [42]:
temp_train = df.iloc[0:len(train_index)]
temp_train = pd.concat((temp_train, target), axis=1)
In [43]:
sns.pairplot(data=temp_train,
y_vars=['SalePrice'],
x_vars=['TotalBsmtSF', '1stFlrSF', '2ndFlrSF', 'BsmtFullBath'])
Out[43]:
In [44]:
X = temp_train[['TotalBsmtSF','2ndFlrSF', 'BsmtFullBath', '1stFlrSF']]
calc_vif(X)
Out[44]:
In [45]:
sns.pairplot(data=temp_train,
y_vars=['SalePrice'],
x_vars=['Total_Bathrooms', 'FullBath', 'HalfBath', 'BsmtHalfBath'])
Out[45]:
In [46]:
X = temp_train[['Total_Bathrooms', 'FullBath', 'HalfBath', 'BsmtHalfBath']]
calc_vif(X)
Out[46]:
In [47]:
sns.pairplot(data=temp_train,
y_vars=['SalePrice'],
x_vars=['OpenPorchSF', '3SsnPorch', 'EnclosedPorch', 'ScreenPorch', 'WoodDeckSF'])
Out[47]:
In [48]:
X = temp_train[['OpenPorchSF', '3SsnPorch', 'EnclosedPorch', 'ScreenPorch','WoodDeckSF']]
calc_vif(X)
Out[48]:
In [49]:
fig, ax =plt.subplots(2,2, figsize=(10,10))
graph = sns.countplot(temp_train['PoolQC'], ax=ax[0,0])
annote_graph(graph)
graph = sns.countplot(temp_train['Utilities'], ax=ax[0,1])
annote_graph(graph)
graph = sns.countplot(temp_train['Street'], ax=ax[1,0])
annote_graph(graph)
fig.show()
fig.delaxes(ax[1][1])
In [50]:
cols_to_remove = ['PoolQC', 'Utilities', 'Street', 'FullBath', '1stFlrSF']
In [51]:
df.drop(cols_to_remove, axis=1, inplace=True)
In [52]:
df.head()
Out[52]:
In [53]:
df.info()
Target variable¶
In [54]:
sns.distplot(target, color='blue')
Out[54]:
In [55]:
from scipy import stats
import numpy as np
z = np.abs(stats.zscore(target))
print(z)
In [56]:
threshold = 3
outlier_index = np.where(z > 3)[0]
print("Number of outliers based on z-score:", len(outlier_index))
The loss function for the regression model is mean_absolute_error which is somewhat robust to outliers compared to mean_squared error which penalizes even a small error which leads to over-estimation of how bad the model is¶
1. Feature Selection using Predictive Power Score¶
In [57]:
import ppscore as pps
In [58]:
temp_train = df.iloc[0:len(train_index)]
temp_train = pd.concat((temp_train, target), axis=1)
temp_train.head()
Out[58]:
In [59]:
col = temp_train.columns#.score will be a column in the matrix below
train_pps = temp_train[col]
pps_mat = pps.matrix(train_pps)
In [60]:
pps_mat[pps_mat['x'] == 'SalePrice'].sort_values(by='ppscore', ascending=False)
Out[60]:
In [61]:
# Selecting features only with ppscore > 0
pps_features = pps_mat.loc[(pps_mat['x'] == 'SalePrice') & (pps_mat['ppscore'] > 0)]['y'].values
In [62]:
print("Number of features using PPS: ",len(pps_features)), pps_features
Out[62]:
Recreating training and testing datasets¶
In [63]:
def create_train_test(df, features, trainIndex):
df = df[features[:-1]]
categorical_cols = df.dtypes[ df.dtypes == 'object' ].index.tolist()
numerical_cols = df.dtypes[ df.dtypes != 'object' ].index.tolist()
cat_dummies = pd.get_dummies(df[categorical_cols], drop_first=True)
df = df.drop(categorical_cols, axis=1)
df = pd.concat([df, cat_dummies], axis=1)
train = df.iloc[0:len(trainIndex)]
test = df.iloc[len(trainIndex):]
print('Shape of Train Dataset:',train.shape)
print('Shape of Test Dataset:',test.shape)
return(train,test)
In [64]:
# Train and Test data with all the features
train_all, test_all = create_train_test(df, df.columns, train_index)
In [65]:
# Train and Test data with features selected by predictive power score
train_pps, test_pps = create_train_test(df, pps_features, train_index)
In [66]:
from sklearn.preprocessing import StandardScaler
def scaler_function(train_data, target):
scaler = StandardScaler()
scaler.fit(train_data)
train_std = scaler.transform(train_data)
train = pd.DataFrame(train_std, columns=train_data.columns)
train = pd.concat((train,target),axis=1)
return scaler, train
In [67]:
scaler_all, train_scaled_all = scaler_function(train_all, target)
scaler_pps, train_scaled_pps = scaler_function(train_pps, target)
In [68]:
train_scaled_pps.shape, train_scaled_all.shape
Out[68]:
Hyperopt¶
In [69]:
from hpsklearn import HyperoptEstimator
from hpsklearn import xgboost_regression, gradient_boosting_regression, extra_trees_regression, ada_boost_regression
from hyperopt import tpe
In [70]:
model_list = [xgboost_regression, gradient_boosting_regression, extra_trees_regression, ada_boost_regression]
model_names = ['xgBoost', 'gradBoost', 'extraTree', 'adaBoost']
In [71]:
from sklearn import metrics
def model_eval(model, y_true, y_pred):
mae = metrics.mean_absolute_error(y_true.values, y_pred)
mse = metrics.mean_squared_error(y_true.values, y_pred)
rmse = np.sqrt(metrics.mean_squared_error(y_true.values, y_pred))
r2score = metrics.r2_score(y_true.values, y_pred)
scores = [mse,rmse, mae,r2score]
return([round(x,5) for x in scores])
In [72]:
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error
import pickle
import os
from pathlib import Path
def modeling(df, target_col, model_list, model_names, test_size=0.3, folder='model_dir', save_models=False):
folder = 'Models/'+folder
p = Path(folder)
p.mkdir(parents=True, exist_ok=True)
# evaluating scores -> r2, rmse, mae
train_scores = []
valid_scores = []
features = df.columns.tolist()
features.remove(target_col)
X = df[features]
y = df[target_col]
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=RANDOM_STATE, test_size=test_size)
for algo, name in zip(model_list,model_names):
model = HyperoptEstimator(regressor=algo(name),
preprocessing=[],
max_evals=N_ITER,
loss_fn=mean_absolute_error,
#algo=tpe.suggest,
trial_timeout=60,
seed=RANDOM_STATE)
model.fit(X_train.values, y_train.values)
#--------- Train Data Evaluation ---------
predictions = model.predict(X_train.values)
scores = model_eval(model= model, y_true= y_train, y_pred= predictions)
train_scores.append(scores)
#--------- Test Data Evaluation ---------
predictions = model.predict(X_test.values)
scores = model_eval(model= model, y_true= y_test, y_pred= predictions)
valid_scores.append(scores)
if save_models:
#print(model.best_model())
pkl_filename = name+'.pkl'
with open(p/pkl_filename, 'wb') as file:
pickle.dump(model, file)
#joblib.dump(model, filename)
train_score_df = pd.DataFrame(train_scores, columns=['mse','rmse', 'mae','r2'], index= model_names)
valid_score_df = pd.DataFrame(valid_scores, columns=['mse','rmse', 'mae','r2'], index= model_names)
return (train_score_df, valid_score_df)
In [74]:
RANDOM_STATE = 42
N_ITER = 40
In [75]:
tr_df, val_df = modeling(df = train_scaled_all,
target_col = target_col,
model_list = model_list,
model_names = model_names,
folder = 'fullData'
)
In [76]:
tr_df
Out[76]:
In [77]:
val_df
Out[77]:
In [78]:
RANDOM_STATE = 42
N_ITER = 40
tr_df_pps, val_df_pps = modeling(df=train_scaled_pps, target_col=target_col,
model_list=model_list, model_names=model_names, folder='pps')
In [79]:
tr_df_pps
Out[79]:
In [80]:
val_df_pps
Out[80]:
In [246]:
# Loading the best performing model -> xgBoost
In [81]:
with open('./Models/fullData/xgBoost.pkl', 'rb') as f:
xgBoost = pickle.load(f)
In [82]:
xgBoost.best_model()
Out[82]:
Final Thoughts¶
In [83]:
n_feat = housePrices.shape[1]-1
n_feat_pps = len(pps_features)-1
print("Number of Independent Features in ORIGINAL dataset", n_feat)
print("Number of Independent Features selected by Predictive Power Score", n_feat_pps)
print("Percentage of features omitted using PPS: ", (n_feat-n_feat_pps)/n_feat * 100)
In case of an regression, the ppscore uses the mean absolute error (MAE) as the underlying evaluation metric (MAE_model).
The best possible score of the MAE is 0 and higher is worse.
As a baseline score, we calculate the MAE of a naive model (MAE_naive) that always predicts the median of the target column.
The PPS is the result of the following normalization (and never smaller than 0):
- PPS = 1 - (MAE_model / MAE_naive)
- Medium article on PPS: Link
- Metrics for the xgBoost model with all features on validation set :
- MAE = 15884.68449 - R-squared = 0.90460 Metrics for the xgBoost model with features selected by PPS :
- MAE = 17013.10227 - R-squared = 0.90928The performance of the model is similar and we have elimated 67.5% of features using PPS.
Test data - Evaluation and Submission¶
In [84]:
test_all.head(3)
Out[84]:
In [85]:
def test_submissions(model, test_data, scaler, model_name, save=False):
test_cols = test_data.columns
test_std = scaler.transform(test_data)
predictions = model.predict(test_std)
d = {'Id': test_index, 'SalePrice': np.round(predictions,3)}
submission = pd.DataFrame(data=d)
if save:
submission.to_csv(('submission_'+model_name+'.csv'),index=False)
return submission
In [92]:
with open('./Models/fullData/xgBoost.pkl', 'rb') as f:
xgBoost_all = pickle.load(f)
In [93]:
xgBoost_all.best_model()
Out[93]:
In [94]:
test_all_predictions = test_submissions(xgBoost_all, test_all, scaler_all, 'xgBoost_test', False)
In [87]:
sns.boxplot(test_all_predictions, showmeans=True)
Out[87]:
In [95]:
test_all_predictions.to_csv('submission_all.csv', index=False)
api.competition_submit('submission_all.csv','All features','house-prices-advanced-regression-techniques')
Out[95]:
In [96]:
with open('./Models/pps/xgBoost.pkl', 'rb') as f:
xgBoost_pps = pickle.load(f)
In [97]:
xgBoost_pps.best_model()
Out[97]:
In [98]:
test_pps_predictions = test_submissions(xgBoost_pps, test_pps, scaler_pps, 'xgBoost_pps', False)
In [99]:
sns.boxplot(test_pps_predictions, showmeans=True)
Out[99]:
In [100]:
test_pps_predictions.to_csv('submission_pps.csv', index=False)
api.competition_submit('submission_pps.csv','PPS features','house-prices-advanced-regression-techniques')
Out[100]: